import gmpy2
from gmpy2 import mpz
import sys

# Miller-Rabin primality test (probabilistic)
def miller_rabin(n, bases):
    if n < 2:
        return False
    if n in (2, 3):
        return True
    if n % 2 == 0:
        return False

    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1

    for a in bases:
        if a % n == 0:
            continue
        x = gmpy2.powmod(a, d, n)
        if x == 1 or x == n - 1:
            continue
        passed = False
        for _ in range(s - 1):
            x = gmpy2.powmod(x, 2, n)
            if x == n - 1:
                passed = True
                break
        if not passed:
            return False
    return True

# Deterministic primality check
def is_prime_deterministic(n):
    return gmpy2.is_prime(n, 0)  # 0 = use deterministic method

# Simple trial division factors (for small pseudoprimes)
def prime_factors(n):
    n = mpz(n)
    factors = []
    while n % 2 == 0:
        factors.append(2)
        n //= 2
    f = 3
    while f * f <= n:
        while n % f == 0:
            factors.append(f)
            n //= f
        f += 2
    if n > 1:
        factors.append(n)
    return factors

# Progress bar
def print_progress(i, total, found):
    bar_length = 40
    progress = i / total
    filled = int(bar_length * progress)
    bar = "#" * filled + "-" * (bar_length - filled)
    percent = int(progress * 100)
    sys.stdout.write(f"\r[{bar}] {percent:3d}% | Found {found}")
    sys.stdout.flush()

if __name__ == "__main__":
    bases = [2, 3]
    limit = 50000000
    print(f"Searching pseudoprimes up to {limit} with bases {bases}...\n")

    count = 0
    total = (limit - 3) // 2 + 1

    for i, n in enumerate(range(3, limit, 2), start=1):
        if i % 1000 == 0 or n == limit:
            print_progress(i, total, count)

        # Step 1: true primality check (deterministic)
        if is_prime_deterministic(n):
            continue  # skip primes

        # Step 2: check if MR test incorrectly says "prime"
        if miller_rabin(n, bases):
            print(f"  Found pseudoprime: {n} = " + " x ".join(map(str, prime_factors(n))))
            count += 1

    sys.stdout.write("\r" + " " * 60 + "\r")  # clear line
    sys.stdout.flush()
    print(f"Search complete. Found {count} pseudoprimes.")